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The Monte Carlo calculation of Renyi entanglement entropies S n of interacting fermions suffers 
from a well-known signal-to-noise problem, even for a large number of situations in which the in¬ 
famous sign problem is absent. A few methods have been proposed to overcome this issue, such 
as ensemble switching and the use of auxiliary partition-function ratios. Here, we present an ap¬ 
proach that builds on the recently proposed free-fermion decomposition method; it incorporates 
entanglement in the probability measure in a natural way; it takes advantage of the hybrid Monte 
Carlo algorithm (an essential tool in lattice quantum chromodynamics and other gauge theories 
with dynamical fermions); and it does not suffer from noise problems. This method displays no 
sign problem for the same cases as other approaches and is therefore useful for a wide variety of 
systems. As a proof of principle, we calculate S 2 for the one-dimensional, half-filled Hubbard model 
and compare with results from exact diagonalization and the free-fermion decomposition method. 
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I. INTRODUCTION 

Topological and quantum-information aspects of con¬ 
densed matter physics, broadly defined to include all 
few- and many-body quantum systems, continue to gain 
increasing attention from a variety of angles, with the 
quantum-mechanical notion of entanglement playing a 
central role [1]. Topological quantum phase transitions, 
in particular, have been shown to bear a direct quanti¬ 
tative connection to the so-called entanglement entropy 
in both its Renyi and von Neumann forms, the entan¬ 
glement spectrum, and other information-related quan¬ 
tities [2]. Thus, the computation of Renyi entangle¬ 
ment entropies S n is currently of vital importance to 
many fields (see e.g. [3-5]), and the challenge of doing 
so in interacting systems, particularly in strongly cou¬ 
pled regimes, must be met. 

To this end, a variety of Monte Carlo (MC) methods 
have recently been put forward to calculate S n efficiently 
(see e.g. [6-13]). As we explain below, one of the cru¬ 
cial steps common to many of the underlying formalisms 
is the so-called replica trick [4], which results in an ex¬ 
pression for S n that consists of a ratio of two partition 
functions. Generally speaking, partition functions them¬ 
selves are challenging objects to compute from the nu¬ 
merical standpoint, as they typically involve terms that 
vary on vastly different numerical scales. In the context 
of stochastic calculations of S n , it is now well understood 
that this complication manifests itself as a signal-to-noise 
problem: the direct estimation of the partition functions, 
followed by the calculation of their ratio, leads to statisti¬ 
cal uncertainties that grow exponentially with the volume 
of the (sub-)system considered (see e.g. [9, 11, 12] for an 
explanation). 


In this work, we present an alternative lattice MC ap¬ 
proach for the calculation of S n . We use a specific case of 
one-dimensional spin-1/2 fermions governed by the Hub¬ 
bard Hamiltonian as an example, which allows us to com¬ 
pare our results with the exact solution as well as with 
other MC methods, but the technique can be generalized 
to arbitrary systems, including those with gauge fields 
and Fermi-Bose mixtures (as long as the so-called sign 
problem is absent, as in any other MC calculation; see 
e.g. [14]). To highlight the generality of our technique, we 
carry out our calculations using the hybrid Monte Carlo 
algorithm (HMC) [15] (see [16] for basic introductions to 
HMC), which is the workhorse of lattice QCD, it is es¬ 
sential in non-perturbative studies of gauge theories with 
dynamical fermions and, more recently, has been used in 
a variety of graphene studies [17 -20]. 


II. FORMALISM 

For the following discussion, we put the system on a 
d-dimensional spatial lattice of side N x . Because we are 
considering such a finite lattice, the single-particle space 
is of finite size N~. We then follow closely the formalism 
of Ref. [9]. 

The n-th Renyi entanglement entropy S n of a sub¬ 
system A of a given quantum system is defined by 

S n = (jT~ lntr (PA), (!) 

where p A is the reduced density matrix of sub-system A 
(i.e. the degrees of freedom of the rest of the system are 
traced over). More concretely, for a system with density 
matrix p, the reduced density matrix is defined via a 
partial trace over the Hilbert space corresponding to the 
complement A of our sub-system as 
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In Ref. [9], Grover derived an auxiliary-field path-integral 
form for p A , from which he showed that S n can be com¬ 
puted using MC methods for a wide variety of systems. 
We summarize those derivations next. In auxiliary- 
field Monte Carlo methods one introduces a Hubbard- 
Stratonovich field a that decouples the fermion species, 
such that the usual density matrix p can be written as a 
path integral: 


where the field integration measure, given by 

z>w= n ^ (s) 

*;=i 

is over the n “replicas” a k of the Hubbard-Stratonovich 
auxiliary field. For convenience, we have included the 
normalization 


e -0H r 

P= —g— = J T>aP[a\p[a], (3) 

for some normalized probability measure P[cr\ deter¬ 
mined by the details of the underlying Hamiltonian (for 
more detail, see [14]). Here, Z = tr[e _/3H ] is the par¬ 
tition function, and p[a] is an auxiliary density matrix 
corresponding to noninteracting particles in the external 
field a. In analogy with this, Grover proved that one may 
decompose the reduced density matrix as 

p A = J V<jP[ct}pa[o-}, (4) 

where 

p A [cr\ = C A [a] exp I -^ctpn^V] “ 

\ 

( 5 ) 

and 

G^fcr] = det(U - G A [a}). (6) 

Here, we have used the restricted Green’s function 
[G^Jo-]]^ , which corresponds to a noninteracting single¬ 
particle Green’s function G(i,j) in the background field 
<7 but such that the arguments i, j only take values in the 
region A (see Ref. [9] and also Ref. [21], where expressions 
for the reduced density matrix of noninteracting systems, 
based on reduced Green’s functions, were first derived). 

Using the above judicious choice of p A [a] , Grover ver¬ 
ified that the expectation value of Cjc\ in the auxiliary 
noninteracting system reproduces the restricted single¬ 
particle Green’s functions, as required. Therefore expec¬ 
tation values of observables supported in the region A 
are reproduced as well. This validates the expression on 
the right-hand side of Eq. (4). 

Using that expression, taking powers of p A results in 
the appearance of multiple auxiliary fields, which we will 
denote below collectively as {<r}. Seemingly an explicit 
manifestation of the replica trick [4], this approach allows 
the trace of the n-th power of p A in Eq. (1) to be recast 
as a multiple field integral over a product of fermion de¬ 
terminants that depend collectively on all the {cr}. In¬ 
deed, for a system of 2iV-component fermions, using a 
Hubbard-Stratonovich transformation that decouples all 
of them, we obtain 

exp ((1 - n)S n ) = J V{cr}P[{a}] Q[{c r}], (7) 



.. 2 N 

Z= Va n det U m [a] (9) 

** m—1 

in the integration measure. The naive probability mea¬ 
sure, given by 


n 2 N 

p iM]=n n detu mWki ( io ) 

fc=1 m= 1 

factorizes entirely across replicas, and is therefore blind 
to entanglement. This factorization is the main reason 
why using P[{cr}] as a MC probability leads to (seem¬ 
ingly) insurmountable signal-to-noise issues, as shown in 
Ref. [9]; it is also why we call it naive (although that 
is by no means our judgement of Ref. [9]). In Eq. (10), 
[o'] is a matrix which encodes the dynamics of the 
?n-th component in the system, namely the kinetic en¬ 
ergy and the form of the interaction after a Hubbard- 
Stratonovich transformation; it also encodes the form of 
the trial state I'F) in ground-state approaches (see e.g. 
Ref. [14]), as is the case in this work. We will take |\H) 
to be a Slater determinant. In finite-temperature ap¬ 
proaches, t/ m [o - ] is obtained by evolving a complete set 
of single-particle states in imaginary time. 

The quantity that contains the essential contributions 
to entanglement is 


2N 

QIW}} = II detM m[ML (ii) 

m—1 


where 




(l G At m[ a k}) X 

fc=1 


n 


i+n 

k =1 


^A,m[ a k} 

H ~ *U4,m[Ofc] 


( 12 ) 


In the above equation, we have used G Am [a k \, which is 
a restricted Green’s function, as previously defined, but 
where we now indicate the fermion component m and 
replica field index k. 

The product Q[{o'}] was identified as playing the role 
of an observable in Ref. [9], which is a natural interpre¬ 
tation in light of Eq. (7), but which we will interpret 
differently below. Note that, for n = 2, no matrix inver¬ 
sion is required in the calculation of Q[{cr}]; for higher 
n, however, there is no obvious way to avoid the inver¬ 
sion of 1 — G Am [a k \. In turn, this would require some 
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kind of numerical regularization technique (see Ref. [13]) 
to avoid the singularities in G A TO [<j fc ], whose eigenvalues 
can be very close to 0 and 1. 

In ground-state approaches, the size of U m [a\ is given 
by the number of particles of the m-th species present in 
the system. In finite-temperature approaches, the size is 
that of the whole single-particle Hilbert space (i.e. the 
size of the lattice). The size of G A m [cr fc ], on the other 
hand, is always given by the number of lattice sites en¬ 
closed by the region A. Note that, separating a factor of 
Z n in the denominator of Eq. (7), an explicit form can be 
identified in the numerator as the result of the so-called 
replica trick [4] (namely a partition function for n copies 
of the system, “glued” together in the region A). 


III. OUR PROPOSED METHOD 

In analogy to the right-hand side of Eq. (7), we in¬ 
troduce an auxiliary parameter 0 < A < 1 and define a 
function T(X;g) via 

T(A; g) = J V{a}P[{a}] Q[{a}; A], (13) 

where we have replaced the dependence of (2[{cr}] on the 
coupling g by setting 

9 -> A 2 ^, (14) 

which defines Q[{cr};A]. Using this definition, it follows 
immediately that T(A; g) satisfies two important con¬ 
straints with physical significance: For A = 0, we find 


Using Eq. (13), 

d In T /' 

— =y Z>MP[W;A]Q[M;A] 

where 

p[W;A] = fpb) p[W] Q K a > ; A], 

and 


2N 


<2[M; A] = tr 


m =1 


MZ\[ {a}} dMm ’ x[{a}] 


m ,A L 


d\ 


(17) 


(18) 


(19) 


Crucially, the dependence on the parameter A enters only 
through the matrix M n . and it is in this way that we pro¬ 
pose to include the entanglement-related correlations in 
the sampling procedure, which is to be carried out using 
P[{cr}; A] as a probability measure. When an even num¬ 
ber 2 N of flavors is considered, and the interactions are 
attractive, det 2Ar U[a\ and Q[{cr}, A] are real and posi¬ 
tive semidefinite for all cr, which means that there is no 
sign problem and P[{cr}; X] above is indeed a well-defined, 
normalized probability measure. 

More concretely, our proposal to calculate S n is to take 
the noninteracting A = 0 point as a reference and com¬ 
pute S n using 

s n = s i 0) + -!— f dX <Q[M; A]), (20) 

1 - n Jo 


where 


--lnr(0 ;g) = S<°\ (15) 

1 — n 

where S„°^ corresponds to a noninteracting system. In¬ 
deed, at A = 0 the quantity Q[{cr};A] does not depend 
on {cr} and factors out of the integral. Thus, regardless 
of the value of g, the function r(0;<?) corresponds to the 
Renyi entropy of a noninteracting system, which can be 
trivially computed with the present formalism. Indeed, 
there is no path integral when interactions are turned 
off, such that the noninteracting result can be computed 
with a single Monte Carlo sample using the formalism 
by Grover mentioned above. It is worth noting at this 
point that the Renyi entropy of a noninteracting system 
has received substantial attention in the last few years. 
Much is known about this quantity for a variety of sys¬ 
tems, in particular in connection with area laws and their 
violation [22]. 

For A = 1, on the other hand, T(A; g) yields the entan¬ 
glement entropy of the fully interacting system: 

—lnT(l;g) = S n . (16) 

1 — n 

Thus, both of these reference points are physically mean¬ 
ingful, one of them is comparatively trivial to obtain, and 
obtaining the other one is our objective. 


{X) = J V{a}P[{a};X}X[{a}}. (21) 

In other words, we obtain an integral form of the in¬ 
teracting Renyi entropy that can be computed using 
any MC method, in particular HMC [15]. The lat¬ 
ter combines molecular dynamics (MD) of the auxiliary 
fields (defining a fictitious auxiliary conjugate momen¬ 
tum) with the Metropolis algorithm, and thus enables si¬ 
multaneous global updates of all the auxiliary fields {cr}. 
As it well known, HMC is a highly efficient sampling 
strategy, particularly when gauge fields are involved (see 
e.g. [14, 15]). The integration of the MD equations of 
motion requires the calculation of the MD force, which 
is given by the functional derivative of P[{cr};A] with 
respect to {cr}, which can be calculated from Eq. (18). 

Our proposal is akin to the so-called coupling-constant 
integration approach of many-body physics, but it differs 
in that we have strategically introduced the A dependence 
only in the entanglement-sensitive determinant Q[{cr}; A] 
ofEq. (11). 

Equation (20) is our main result and defines our 
method. An essential point is that the expectation that 
appears above is taken with respect to the probabil¬ 
ity measure P[{cr};A] , which incorporates the correla¬ 
tions that account for entanglement. In stark contrast 
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to the naive MC probability P[{er}], which corresponds 
to statistically independent copies of the Hubbard- 
Stratonovich field, this distribution does not display the 
decoupling responsible for the signal-to-noise problem 
mentioned above. 

In practice, using Eq. (20) requires MC calculations 
to evaluate (Q[{cr};A]) as a function of A, followed by 
numerical integration over A. We find that (Q[{a};A]) 
is a smooth function of A that vanishes at A = 0 and 
presents most of its features close to A = 1 (further details 
below). We therefore carry out the numerical integration 
using the Gauss-Legendre quadrature method [23]. It 
should also be pointed out that the stochastic evaluation 
of (Q[{ct};A]), for fixed subregion A, can be expected to 
feature roughly symmetric fluctuations about the mean. 
Therefore, after integrating over A, the statistical effects 
on the entropy are reduced (as we show empirically in 
the Results section). 

A few remarks are in order regarding the auxiliary pa¬ 
rameter A. First, we could have performed the replace¬ 
ment g —> A 2 g everywhere, i.e. not only in Q but also in P 
(and its normalization Z). This would have led to three 
terms in the derivative of PQ with respect to A, two of 
which would come from P (recall P is normalized) and 
feature different signs and a rather indirect connection 
to S n (recall P factorizes across replicas). Our approach 
avoids this extra complication by focusing on the entan¬ 
glement part of the path integral (i.e. Q). Second, we 
could have used \ x instead of A 2 , where x does not have 
to be an integer (although it would be rather inconve¬ 
nient to make it less than 1). This is indeed a possibility 
and it allows for further optimization than pursued here. 
In the remainder of this work we set x = 2, as above. 
Finally, the required calculations for different values of 
A are completely independent from one another and can 
therefore be performed in parallel with essentially perfect 
scaling (up to the final data gathering and quadrature). 


IV. RELATION TO OTHER METHODS 


one introduces auxiliary ratios of the partition function 
corresponding to systems whose configuration spaces are 
only marginally dissimilar. In other words, one writes 


exp((l-n)SJ 


2- A,n 

Z n 


-2 A,n -2 A-l,n 

•2)4-1,n -2)4-2,n 



-^l , n 

Z n ’ 
(23) 


where the auxiliary ratios Z A _ in /Z A _ i _ ln are chosen 
to correspond to subsystems of similar size and shape 
(e.g. such that their linear dimension differs by one lat¬ 
tice point). In this way, each of the auxiliary ratios can 
be expected to not differ significantly from unity. With 
enough intermediate ratios, calculations can be carried 
out in a stable fashion at the price of calculating a po¬ 
tentially large number of ratios. 

In the method we propose here, the parameter A plays 
the role of the varying region size A of the ratio trick. 
Indeed, using Eq. (20), we may schematically write 


exp ((1 — n)(S n - Sl 0) )) = exp (AA (Q[M; A])) , 

A=0 

. (24) 

where any discretization AA of the exponent inside 
the product yields a telescopic sequence of ratios as in 
Eq. (23). As long as (Q[{c};A[) is regular in A, which 
we find to be the case, our auxiliary factors can be made 
to be arbitrarily close to unity at the cost of (at most) 
linear scaling in computation time. 


V. RESULTS 

We test our algorithm by computing the second Renyi 
entropy S 2 for one-dimensional, ten-site, half-filled Hub¬ 
bard models with periodic boundary conditions. The 
Hamiltonian we used is 

H = ~t ( 4 , A* + cj ;S c iiS ) + U ^n it n 4 , (25) 

a,(ij) i 


Our approach is very similar to the temperature- 
integration method of Ref. [6], but is closer in nature 
to the ratio trick (and similar) of Refs. [7, 12]. As above, 
the calculation starts from the replica trick of Calabrese 
and Cardy [4], i.e. 

exp((l - n)S n ) = -2^, (22) 

where Z A is the partition function of n copies of the 
system “glued” together in the region A. Typically, Z A n 
and Z n can be very different from each other in magni¬ 
tude, particularly if S n is large (as is typically the case 
for large sub-system sizes). Therefore, computing the 
above partition functions separately (and stochastically) 
and then attempting to evaluate the ratio is likely to 
yield a large statistical uncertainty. A way around this 
problem is to use the ratio (or increment) trick, whereby 


where the first sum includes two fermion flavors s =f, 4- 
and nearest-neighbor pairs. To carry out our tests, we 
implemented a symmetric Trotter-Suzuki decomposition 
of the Boltzmann weight, with an imaginary-time dis¬ 
cretization parameter r = 0.05 (in lattice units). The 
full length of the imaginary-time direction was at most 
/3 = 5 (i.e. we used 100 imaginary-time lattice points). 
The interaction factor in the Trotter-Suzuki decomposi¬ 
tion was addressed, as anticipated in a previous section, 
by introducing a replica auxiliary field a for each power 
of the reduced density matrix. This insertion was ac¬ 
complished via a Hubbard-Stratonovich transformation, 
which we chose to be of a continuous but compact form 
(see Ref. [14]). 

Figure 1 plots (Q[{cr}; A]) as a function of both A and 
the subregion size La for four values of the coupling. 
We note that surfaces corresponding to weak couplings 
show much less fluctuation in both parameters than their 
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strong-coupling counterparts. This uniformity implies detailed above, 
that for weakly coupled systems, even at large La, a 
coarse A discretization may yield good estimates of S n . 

Conversely, strongly coupled systems are, not unexpect¬ 
edly, more computationally demanding. Comparison with exact diagonalization results 

and a first look at statistical effects 


-dS 2 /dA 


U/t = 0.5 
U/t = 1.0 
U/t = 2.0 
U/t = 4.0 



0.4 
0.6L a /L 


FIG. 1: (color online) ffybrid Monte Carlo results for the 
source with U/t = 0.5,1.0, 2.0 and 4.0 (from bottom to top 
at A = 1) as functions of the parameter A and the region size 
La, for N x = 10 sites. 

For small subsystems, in addition to relatively small 
variation, the majority of the deviation from the nonin¬ 
teracting entropy is accumulated at large A and appears 
with mostly uniform signature. Much to the contrary, for 
larger subsystems, intermediate values of A correspond 
to a region of parameter space that contributes opposite- 
sign corrections to the entropy, which yields increasing 
uncertainty as a function of the subsystem size. This 
effect is most clearly seen for the curve corresponding 
to the case where the subregion constitutes the entire 
system. In that case, S 2 is zero regardless of the cou¬ 
pling, which means that the integral over A must van¬ 
ish. This happens by a precise cancelation that must be 
captured by the numerical integration procedure. Given 
that the features of (Q[{a}; A]) are roughly concentrated 
in 0.5 < A < 1, we chose the Gauss-Legendre quadrature 
method to carry out the integral in a precise fashion. Us¬ 
ing N x = 20 points in the interval [0,1] (i.e. 40 points 
in the defining interval [—1,1] using an even extension 
of the integrand), we find that, for the parameter values 
studied here, the systematic effects associated with A are 
smaller than the statistical uncertainty. 

Our experience, as detailed above, indicates that the 
features of (Q[{<r};A]) are generic: they vary in ampli¬ 
tude with the coupling but are largely insensitive to the 
overall system size, and generally behave in a benign way 
as a function of A and the sub-system size. Therefore the 
A integration does not contribute to the scaling of the 
computation time vs. system size beyond a prefactor. 
Next, we present our results upon integrating over A as 


In Fig. 2, we show results for a system of size L = N x i , 
where N x = 10 sites and l = 1 is the lattice spacing (as in 
conventional Hubbard-model studies). Our results cover 
couplings U/t = 0.5, 1, 2, and 4, and subsystems of size 
L a = 1 — 10, which we compare to the results of Ref. [9]. 
The agreement is quite satisfactory. 



FIG. 2: (color online) Results for the ten-site Hubbard model 
with U/t = 0.5, 1.0, 2.0, and 4.0 (solid lines from top to bot¬ 
tom). The results for the noninteracting case are shown with 
a dashed line. Hybrid Monte Carlo answers with numerical 
uncertainties for 7,500 decorrelated samples, shown with er¬ 
ror bars. Exact-diagonalization results of Ref. [9] shown with 
lines, except for the U/t = 4.0 case, where the lines join the 
central values of our results and are provided to guide the eye. 


To better understand the size of the statistical effects, 
we also show how our results vary with the number of MC 
samples in Fig. 3. In Fig. 3 we see that the 25,000 sam¬ 
ples collected were well beyond what was needed: half as 
many would have already given excellent results. These 
results show that, by including entanglement-sensitive 
contributions into the probability measure, our approach 
circumvents the signal-to-noise problem mentioned in the 
introduction. Below we elaborate more explicitly on sta¬ 
tistical effects and said problem, and show concrete nu¬ 
merical examples of how it arises in practice. 

In Fig. 4, we show the overall statistical uncertainty 
AS 2 in our estimates of S 2 as a function of the number 
of samples N s , for the coupling strengths and subsystem 
sizes studied above. A S 2 was computed by using the 
envelope determined by the MC statistical uncertainties 
in (Q[{v}', A]) as a function of A. While A S 2 grows with 

1 /2 

the sub-system size, its N s ' scaling remains constant as 
the number of samples is increased. 
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FIG. 3: (color online) Second Renyi entropy (scaled to the 
noninteracting result) as a function of the number of samples 
N g for coupling U/t = 0.5, 2.0, and 4.0 shown top to bottom. 
Within a few thousand samples, we observe that the results 
have stabilized to within 1-2%. 


B. Comparison with naive free-fermion 
decomposition method 
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FIG. 4: (color online) Relative statistical uncertainty of the 
second Renyi entropy (propagated from the standard devia¬ 
tion in the uncertainties on its A derivative) as a function of 
the number of samples N s for couplings U/t = 0.5, 2.0, and 
4.0 shown top to bottom. The symbols and colors correspond 
to those utilized in Fig. 3. 


Figure 5 again shows our results for the ten-site Hub¬ 
bard model using 7,500 decorrelated samples for each 
value of A, this time compared with the naive free-fermion 
decomposition method using 75,000 samples. Such an in¬ 
creased number of samples for the naive method was cho¬ 
sen to provide a more fair comparison with our method, 
considering that the latter requires a MC calculation for 
each value of A. We used 20 A points but, as explained 
in more detail below, roughly half of the A points require 
only a small number of samples. 

The statistical uncertainties for the naive method do 
not encompass the expected answers for many of the 
points, which is a symptom of an “overlap” problem, 
i.e. the probability measure employed bears little corre¬ 
lation with the observable, as mentioned above (see also 
Ref. [24]). This is the same issue as the signal-to-noise 
problem referred to above. 

To illustrate this situation more precisely, we show in 
Fig. 6 a histogram of <3[{cr}] [see Eq. (7)] for U/t = 2.0 
and La/L = 0.8. Even using a logarithmic vertical scale, 


the distribution displays a long tail that extends across 
multiple orders of magnitude. We find that the distri¬ 
bution is approximately of the log-normal type (i.e. its 
logarithm is approximately distributed as a gaussian, as 
shown in the inset of Fig. 6); this is the challenge faced 
when attempting to determine the average of Q[{cr}] with 
good precision implementing the free-fermion decompo¬ 
sition of Ref. [9] at face value. Moreover, we expect these 
features to worsen in larger systems, higher dimensions, 
and stronger couplings, as the matrices involved become 
more ill-conditioned. 

Notably, it is the logarithm of the expectation value of 
<3[{er}] that determines the entanglement entropy, which 
could then be obtained using the cumulant expansion. 
However, it is a priori entirely unknown whether such an 
expansion would converge, i.e. we do not know to what 
extent this distribution deviates from gaussianity. 

The log-normality referred to above has been associ¬ 
ated with the auxiliary field representation of the inter¬ 
action. In such external fields cr, the orbitals of the trial 
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FIG. 5: (color online) Hybrid Monte Carlo results (solid) for 
the ten-site Hubbard model with U/t = 0.5, 2.0, and 4.0 
with numerical uncertainties for 7,500 decorrelated samples 
(for each value of A) compared with results from the naive 
free-fermion decomposition method (crosses with dashed error 
bars) with 75,000 decorrelated samples. 


wavefunction diffuse much like electrons in a disordered 
medium, and the stronger the interaction (or the lower 
the temperature) the heavier the tail becomes in the dis¬ 
tribution of Q[{er}]. This effect was noticed relatively 
recently in Ref. [24], and it appears to be quite ubiqui¬ 
tous. It was then shown, phenomenologically, that many 
signal-to-noise problems are characterized by the heavy 
tail of a lognormal distribution (see also, Ref. [25]). 


C. Statistical behavior as a function of coupling, 
region size, and auxiliary parameter 

In Fig. 7 we show the statistical distribution of our re¬ 
sults for the A derivative, for several couplings. The dis¬ 
tributions we observe are approximately gaussian (they 
decay faster than linearly in a log scale), except for the 
strongest coupling we studied U/t = 4.0 where, not un¬ 
expectedly, the distribution becomes more asymmetric 
and develops heavier tails relative to its weaker-coupling 



Q 

FIG. 6: (color online) Distribution of the observable Q[{u}] 
of the naive free-fermion decomposition method [implemented 
via Eq. (7)] for U/t = 2.0 and L A /L = 0.8. Note that Q[{u}] 
is a non-negative quantity. The long tail (main plot; note 
logarithmic scale in vertical axis) extends beyond the range 
shown and is approximately a log-normal distribution, i.e. 
InQ[{cr}] is roughly distributed as a gaussian (inset). 


counterparts. 



FIG. 7: (color online) Histogram showing the statistical dis¬ 
tribution of our results for dS 2 /dX for L A /L = 0.8, A ~ 0.83, 
and U/t = 0.5,1.0, 2.0,4.0. The results for different couplings 
have been shifted for display purposes, but the scale is the 
same for each of them. This illustrates that, even though 
our method addresses the original signal-to-noise issue, strong 
couplings remain more challenging that weak couplings. 


In Fig. 8 we plot the statistical distribution of our re¬ 
sults for the A derivative at fixed region size and cou¬ 
pling, but varying A. As claimed above, the chosen 
parametrization requires considerably less MC samples 
at low A than at high A, as the width of the distributions 
is much smaller for the former than for the latter. 

Finally, in Fig. 9 we show the same distribution as 
above, but as a function of subregion size. As expected, 













VI. SUMMARY AND CONCLUSIONS 
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FIG. 8: (color online) Histogram showing the statistical dis¬ 
tribution of our results for dS 2 /d\ for L A /L — 0.8, A ~ 0.225, 
0.369, 0.833, 0.991, and U/t = 2.0. The results for different 
couplings have been shifted for display purposes, but the scale 
is the same for each of them. This illustrates that low values 
of A require less samples than larger ones in order to deter¬ 
mine (<2[{n}; A]) with good precision. 

large subregions are more challenging, but the overall 
shape of the distributions is very well controlled: it is 
close to gaussian in that its tails decay faster than linearly 
in a log scale. 



FIG. 9: (color online) Histogram showing the statistical dis¬ 
tribution of our results for dS 2 /d\ for L A /L =0.1, 0.3, 0.5, 
0.8, A ~ 0.83, and U/t = 2.0. The results for different cou¬ 
plings have been shifted for display purposes, but the scale is 
the same for each of them. 


In this work, we have put forward an alternative MC 
approach to the calculation of the Renyi entanglement 
entropy of many-fermion systems. As an essential fea¬ 
ture of our method, we compute the derivative of the 
entanglement entropy with respect to an auxiliary pa¬ 
rameter and integrate afterwards. We have shown that 
such a derivative can be computed using a MC approach 
without signal-to-noise issues, as the resulting expres¬ 
sion yields a probability measure that does not factor 
across replicas and accounts for entanglement in the MC 
sampling procedure in a natural way. The subsequent 
numerical integration can be carried out efficiently via 
Gauss-Legendre quadrature. 

As a proof of principle, we have presented results for 
S 2 for the ID Hubbard model at half filling for different 
coupling strengths and compared with answers obtained 
by exact diagonalization. Our calculations show that the 
statistical uncertainties are well controlled, as we have 
shown in numerous plots and histograms. Although we 
have not run into numerical stability issues in our tests, 
we anticipate that those may appear in the form dis¬ 
cussed in Ref. [13]. Our approach is just as general as the 
one proposed in Ref. [9]. In particular, it can be straight¬ 
forwardly generalized to finite temperature as well as to 
relativistic systems, in particular those with gauge fields 
such as QED and QCD, or any SU (N) gauge theory. 
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